% Improve this code
clear all
clc, clf
N = 20;
x = linspace(-pi,pi,N);
y = exp(x).*sin(x);
yp1 = -exp(-pi);
ypn = -exp(pi);
% h = diff(x);
h = 2*pi/(N-1);
% lam = [1, ]
A = diag([1/2*ones(1,N-2),1],-1) +diag(2*ones(1,N),0) + diag([1 1/2*ones(1,N-2)],1);
d = [6/h*((y(2)-y(1))/h-yp1), 6/(h+h)*((y(3:N)-y(2:N-1))/h-(y(2:N-1)-y(1:N-2))/h) ,6/h*(ypn-(y(N)-y(N-1))/h)];
M = A/d;
M = M';
alph = y;
gam = M/2;
bet = (y(2:N)-y(1:N-1))/h-(2*M(1:N-1)+M(2:N))/6*h;
del = (M(2:N)-M(1:N-1))/(6*h);

xx = linspace(-pi,pi,100);
yx = exp(xx).*sin(xx);
plot(xx,yx)
hold on

for j=1:N-1
    xx = linspace(x(j),x(j+1),20);
   y1 = alph(j) + bet(j)*(xx-x(j))+gam(j)*(xx-x(j)).^2+del(j)*(xx-x(j)).^3;
   plot(xx,y1)
   hold on
end
figure(2)
for j=1:N-1
    xx = linspace(x(j),x(j+1),20);
    yp = exp(xx).*(sin(xx)+cos(xx));
   y2 =  bet(j)+2*gam(j)*(xx-x(j))+3*del(j)*(xx-x(j)).^2;
   plot(xx,y2,'ro',xx,yp)
   hold on
end